En esta práctica utilizaremos los registros de la especie Leopardus wiedii, la paquetería ENMEval 2.0 (Kass et. al 2025) y el algoritmo MaxEnt (Phillips et al. 2006) para explorar diferentes métodos de partición de datos en la modelación de distribución de especies. El objetivo de esta práctica es comparar tipos de partición de datos, evaluando cómo cada enfoque maneja la autocorrelación espacial entre los datos de entrenamiento y validación, así como su efecto en las métricas.
# Leer datos biológicos
datos <- read.csv("/Users/jimenaburgos/Downloads/practica_8p1/train_lw.csv")
occs <- datos[,c("long","lat")]
# Leer los datos ambientales
setwd("/Users/jimenaburgos/Downloads/CursoNichos_2026-1/Practica2_Evaluation/envs")
envs.files <- list.files(pattern='.tif', full.names=TRUE)
# Hacemos el raster stack
envs <- rast(envs.files)
# Graficamos para ver los registros
plot(envs[[3]])
points(occs, pch=19)
### Cargamos una capa de sesgo para dirigir los puntos de background
biasfile = rast("/Users/jimenaburgos/Downloads/CursoNichos_2026-1/Practica2_Evaluation/biaslayer_HREFBiv.tif")
biaslayer <- biasfile
# Es necesario que no existan NAs en el raster
# Por lo que asignamos un valor de 0 a los NAs
values(biasfile)[is.na(values(biasfile))] <- 0
n.bg <- 1000 # este sera la cantidad de puntos de background
bg <- xyFromCell(biasfile,
sample(ncell(biasfile),
n.bg, prob=values(biasfile))) %>% as.data.frame()
colnames(bg) <- colnames(occs)
head(bg)
## long lat
## 1 -87.02083 20.43750
## 2 -90.43750 19.89583
## 3 -87.81250 18.97917
## 4 -87.89583 20.81250
## 5 -91.06250 18.18750
## 6 -87.39583 21.22917
# Graficamos los puntos de background sobre la capa de sesgo
plot(biaslayer, main = "Capa de sesgo")
points(bg, pch=3 , cex = 0.4, col="white")
# Definimos una función personalizada para implementar la pROC
proc <- function(vars) {
proc <- kuenm::kuenm_proc(vars$occs.val.pred, c(vars$bg.train.pred, vars$bg.val.pred))
out <- data.frame(proc_auc_ratio = proc$pROC_summary[1],
proc_pval = proc$pROC_summary[2], row.names = NULL)
return(out)
}
Los dos métodos siguientes no tienen en cuenta la autocorrelación espacial entre los registros de validación y entrenamiento. Por lo que, es importante considerar que la partición aleatoria puede llevar a la sobreestimación del modelo si los registros de presencia de calibración y evaluación están próximas entre sí, debido a que las localidades utilizadas para evaluar el modelo no son independientes de las utilizadas para calibrarlo.
Principalmente, al trabajar con conjuntos de datos relativamente pequeños (e.g. 25 registros de presencia), los usuarios pueden optar la validación cruzada k-fold, en el que el número de bins (k) o grupos es igual al número de registros de presencia (n) (Pearson et al., 2007; Shcheglovitova y Anderson, 2013). Esto se conoce como partición jackknife o “dejar uno fuera” (Hastie et al., 2009).
jack <- get.jackknife(occs, bg)
evalplot.grps(pts = occs, pts.grp = jack$occs.grp, envs = envs)
# Corremos el modelo
e.mx_jack <- ENMevaluate(occs = occs, #Registros de presencia
envs = envs, #Capas ambientales
bg = bg, #Puntos de background
algorithm = 'maxnet', #Algoritmo
partitions = 'jackknife', #Método de partición
tune.args = list(fc = c("L","LQ","LQP"), rm = 1),
user.eval = proc) #Métrica adicional (pROC)
## | | | 0% | |======================= | 33% | |=============================================== | 67% | |======================================================================| 100%
El método de «k-fold aleatorio» divide aleatoriamente los registros de presencia en un número de (k) bins o grupos especificado por el usuario (Hastie et al., 2009). Este método es ideal utilizarlo con conjuntos de registros de presencia grandes.
rand <- get.randomkfold(occs, bg, k = 5) #En este caso elegimos 5
evalplot.grps(pts = occs, pts.grp = rand$occs.grp, envs = envs)
# Corremos el modelo
e.mx_kfold <- ENMevaluate(occs = occs,
envs = envs,
bg = bg,
algorithm = 'maxnet',
partitions = 'randomkfold',
tune.args = list(fc = c("L","LQ","LQP"), rm = 1),
user.eval = proc)
## | | | 0% | |======================= | 33% | |=============================================== | 67% | |======================================================================| 100%
Los métodos de partición son variaciones de las particiones geográficamente estructuradas (Radosavljevic y Anderson, 2014). Básicamente, estos métodos dividen los registros de presencia y background en grupos de evaluación basados en reglas espaciales. El objetivo es reducir la autocorrelación espacial entre los registros utilizados en la validación y el entrenamiento, evitando inflar el rendimiento del modelo, al menos para los conjuntos de datos resultantes de un muestreo sesgado (Veloz 2009, Wenger y Olden 2012, Roberts et al. 2017).
# El método de "bloque" particiona los datos según la mediana de la latitud
# y longitud que dividen los registros de presencia en cuatro grupos espaciales de
# igual número o lo más próximos posible.
block <- get.block(occs, bg, orientation = "lat_lon")
# El objeto resultante es una lista de dos vectores que proporcionan la
# designación del bloque o bin para cada registro de presencia y background
table(block$occs.grp)
##
## 1 2 3 4
## 10 10 10 9
# La primera dirección divide los registros de presencia en dos grupos, y la
# segunda divide cada uno de ellos en dos grupos adicionales, resultando en cuatro grupos.
# Tanto los registros de presencia como los de background se asignan a cada uno de
# los cuatro bloques o bins, según su posición con respecto a estas líneas
# bloques de registros de presencia
evalplot.grps(pts = occs, pts.grp = block$occs.grp, envs = envs) +
ggplot2::ggtitle("Spatial block partitions: occurrences")
# bloques de background
evalplot.grps(pts = bg, pts.grp = block$bg.grp, envs = envs) +
ggplot2::ggtitle("Spatial block partitions: background")
# Corremos el modelo
e.mx_block <- ENMevaluate(occs = occs,
envs = envs,
bg = bg,
algorithm = 'maxnet',
partitions = 'block',
tune.args = list(fc = c("L","LQ","LQP"), rm = 1),
user.eval = proc)
## | | | 0% | |======================= | 33% | |=============================================== | 67% | |======================================================================| 100%
Estos generan cuadrículas como un tablero de ajedrez a lo largo de la extensión del área de calibración y dividen los registros de presencia en grupos según su ubicación en el tablero. A diferencia del método de bloques, ambos métodos de tablero de ajedrez subdividen el área de calibración equitativamente, aunque no garantizan el mismo número de localidades de ocurrencia en cada bin o grupo.
El método de tablero de ajedrez básico divide los puntos en k = 2 grupos utilizando un patrón de tablero de ajedrez simple.
cb1 <- get.checkerboard(occs, envs, bg,
aggregation.factor=30) #puedes cambiar
#este parámetro para crear "cajas" más grandes.
# visualizamos presencias
evalplot.grps(pts = occs, pts.grp = cb1$occs.grp, envs = envs)
# visualizamos background
evalplot.grps(pts = bg, pts.grp = cb1$bg.grp, envs = envs)
# Corremos el modelo
e.mx_check1 <- ENMevaluate(occs = occs,
envs = envs,
bg = bg,
algorithm = 'maxnet',
tune.args = list(fc = c("L","LQ","LQP"), rm = 1),
partitions = "user",
user.grp = cb1,
user.eval = proc)
## | | | 0% | |======================= | 33% | |=============================================== | 67% | |======================================================================| 100%
El método de tablero de Checkerboard jerárquico genera k = 4 grupos espaciales mediante un enfoque de dos escalas anidadas. Este método aplica dos factores de agregación secuenciales: primero crea una grilla gruesa que divide el área de estudio, y luego subdivide cada celda de esta grilla en celdas más finas, generando un patrón espacial jerárquico.
A diferencia del checkerboard básico que produce solo 2 grupos, este método ofrece una partición más granular con 4 grupos, manteniendo la separación geográfica pero proporcionando mayor flexibilidad para conjuntos de datos grandes o áreas de estudio complejas. Los registros de presencia y background se asignan a cada grupo según su ubicación en el patrón jerárquico resultante.
cb2 <- get.checkerboard(occs, envs, bg, aggregation.factor=c(10,10))
# visualizamos presencias
evalplot.grps(pts = occs, pts.grp = cb2$occs.grp, envs = envs)
# visualizamos background
evalplot.grps(pts = bg, pts.grp = cb2$bg.grp, envs = envs)
# Corremos el modelo
e.mx_check2 <- ENMevaluate(occs = occs,
envs = envs,
bg = bg,
algorithm = 'maxnet',
tune.args = list(fc = c("L","LQ","LQP"), rm = 1),
partitions = "user",
user.grp = cb2,
user.eval = proc)
## | | | 0% | |======================= | 33% | |=============================================== | 67% | |======================================================================| 100%
# Podemos ver los resultados de cualquiera de nuestros modelos
# vemos el de jackknife
e.mx_check1@results
## fc rm tune.args auc.train cbi.train auc.diff.avg auc.diff.sd auc.val.avg
## 1 L 1 fc.L_rm.1 0.6900187 0.771 0.06782160 0.07391116 0.6574025
## 2 LQ 1 fc.LQ_rm.1 0.6902389 0.895 0.06133692 0.06881793 0.6577346
## 3 LQP 1 fc.LQP_rm.1 0.6990200 0.901 0.06518888 0.08016369 0.6558617
## auc.val.sd cbi.val.avg cbi.val.sd or.10p.avg or.10p.sd or.mtp.avg or.mtp.sd
## 1 0.05742105 0.4685 0.01343503 0.2045455 0.2892710 0.11363636 0.16070609
## 2 0.05403649 0.4565 0.05161880 0.1818182 0.2571297 0.09090909 0.12856487
## 3 0.06439568 0.3720 0.20364675 0.2045455 0.2892710 0.06818182 0.09642365
## proc_auc_ratio.avg proc_auc_ratio.sd proc_pval.avg proc_pval.sd AICc
## 1 1.167737 0.07570549 0 0 897.1241
## 2 1.177875 0.05026016 0 0 900.6929
## 3 1.159173 0.08067267 0 0 902.5415
## delta.AICc w.AIC ncoef
## 1 0.000000 0.81002561 4
## 2 3.568711 0.13600752 5
## 3 5.417391 0.05396688 6
# Asignamos los resultados de la evaluacion en un dataframe
res_block <- eval.results(e.mx_block)
res_check1 <- eval.results(e.mx_check1)
res_check2 <- eval.results(e.mx_check2)
res_jack <- eval.results(e.mx_jack)
res_kfold <- eval.results(e.mx_kfold)
# Agregamos una columna para distinguir los tipos de partición
res_block["Tipo_particion"] <- "Bloque"
res_check1["Tipo_particion"] <- "CheckerboardBasico"
res_check2["Tipo_particion"] <- "CheckerboardJerarquico"
res_jack["Tipo_particion"] <- "Jackknife"
res_kfold["Tipo_particion"] <- "Random-k-fold"
# Unimos los dataframes
data_results <- rbind(res_block, res_check1, res_check2, res_jack, res_kfold)
data_long <- data_results %>%
# Seleccionar las columnas necesarias
select(Tipo_particion, or.10p.avg, auc.val.avg ) %>%
# Convertir a formato largo directamente (sin agrupar ni promediar)
pivot_longer(
cols = c(or.10p.avg, auc.val.avg),
names_to = "Metrica",
values_to = "Valor"
) %>%
# Renombrar las métricas para que se vean mejor en la gráfica
mutate(
Metrica = case_when(
Metrica == "or.10p.avg" ~ "OR",
Metrica == "auc.val.avg" ~ "AUC",
TRUE ~ Metrica
)
)
# Visualizamos las métricas de los distintos tipos de partición
ggplot(data_long, aes(x = Tipo_particion, y = Valor, color = Tipo_particion, fill = Tipo_particion)) +
geom_boxplot(alpha = 0.7) +
facet_wrap(~Metrica, scales = "free_y") +
theme_minimal() +
labs(
x = "Tipo de partición",
y = "Valor") +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none")
# Utilizamos un protocolo secuencial para elegir un modelo de cada tipo de partición:
# la menor tasa de omisión y el máximo AUC
opt.seq_block <- res_block %>%
filter(or.10p.avg == min(or.10p.avg)) %>%
filter(auc.val.avg == max(auc.val.avg))
opt.seq_check1 <- res_check1 %>%
filter(or.10p.avg == min(or.10p.avg)) %>%
filter(auc.val.avg== max(auc.val.avg))
opt.seq_check2 <- res_check2 %>%
filter(or.10p.avg == min(or.10p.avg)) %>%
filter(auc.val.avg == max(auc.val.avg))
opt.seq_jack <- res_jack %>%
filter(or.10p.avg == min(or.10p.avg)) %>%
filter(auc.val.avg == max(auc.val.avg))
opt.seq_kfold <- res_kfold %>%
filter(or.10p.avg == min(or.10p.avg)) %>%
filter(auc.val.avg == max(auc.val.avg))
# Ahora seleccionamos solo ese modelo de nuestro conjunto de modelos
mod.seq_block <- eval.models(e.mx_block)[[opt.seq_block$tune.args]]
mod.seq_check1 <- eval.models(e.mx_check1)[[opt.seq_check1$tune.args]]
mod.seq_check2 <- eval.models(e.mx_check2)[[opt.seq_check2$tune.args]]
mod.seq_jack <- eval.models(e.mx_jack)[[opt.seq_jack$tune.args]]
mod.seq_kfold <- eval.models(e.mx_kfold)[[opt.seq_kfold$tune.args]]
# Finalmente seleccionamos una predicción
pred.seq_block <- eval.predictions(e.mx_block)[[as.character(opt.seq_block$tune.args)]]
pred.seq_check1 <- eval.predictions(e.mx_check1)[[as.character(opt.seq_check1$tune.args)]]
pred.seq_check2 <- eval.predictions(e.mx_check2)[[as.character(opt.seq_check2$tune.args)]]
pred.seq_jack <- eval.predictions(e.mx_jack)[[as.character(opt.seq_jack$tune.args)]]
pred.seq_kfold <- eval.predictions(e.mx_kfold)[[as.character(opt.seq_kfold$tune.args)]]
# Graficamos los modelos seleccionados de cada tipo de partición
par(mfrow=c(1, 3))
plot(pred.seq_block, main="Block")
plot(pred.seq_kfold, main="Random-k-fold")
plot(pred.seq_check1, main="Checkerboard básico")
par(mfrow=c(1, 2))
plot(pred.seq_check2, main="Checkerboard jerárquico")
plot(pred.seq_jack, main="Jackknife")